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We extend our previous derivation of an exact expression for the leading-order (LO) gluon distri- 
bution function G{x,Q^) = xg{x,Q^) from the DGLAP evolution equation for the proton structure 
. . . function _F2'''(3^i Q^) for deep inelastic 7*p scattering to include the effects of heavy-quark masses. We 

0^ I derive the equation for G{x,Q^) in two different ways, first using our original differential-equation 

. method, and then using a new method based on Laplace transforms. The results do not require the 

' use of the gluon evolution equation, or, to good approximation, knowledge of the individual quark 

distributions. Given an analytic expression that successfully reproduces the known experimental 
O . data for -F^^(x, Q^) in a domain 'D{x, Q^) — where Xmin{Q^) < x < a;niax(Q^), Qniin < < Qmax of 

^ ' the Bjorken variable x and the virtuality — G{x, Q^) is uniquely determined in the same domain. 

' As an application of the method, we construct a new global parametrization of the complete set 

of ZEUS data on F:^^{x, Q'^), and use this to determine the 5 quark gluon distribution, G{x,Q ), 
for massless u, d, s and massive c, b quarks and discuss the mass effects evident in the result. We 
compare these results to the gluon distributions for CTEQ6L, and in the domain D{x, Q^) where 
they should agree, they do not; the discrepancy is due to the fact that the CTEQ6L results do 
f~| . not give an accurate description of the ZEUS F2^{x,Q^) experimental data. We emphasize that 

^ l' our method for obtaining the the LO gluon distribution connects G{x,Q^) directly to the proton 

structure function without either the need for individual parton distributions or the gluon evolution 
equation. 



PACS numbers: 13.85.Hd,12.38.Bx,12.38.-t,13.60.Hb 
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> ' I. INTRODUCTION 

(N 
. 

, The quark and gluon distributions in hadrons play a key role in our understanding of Standard Model processes, 
' in our predictions for such processes at accelerators, and in our searches for new physics. In particular, accurate 
' knowledge of gluon distribution functions at small Bjorken x will play a vital role in estimating backgrounds, and 
. hence, our ability to search for new physics at the Large Hadron Collider. 

' The gluon and quark distribution functions have traditionally been determined simultaneously by fitting expcri- 
. mental data (mainly at small x) on the proton structure function _F^^(a;, Q^) measured in deep inelastic ep (or j*p) 
^ ' scattering, over a large domain of values of x and . The process starts with an initial Qq, typically in the 1 to 2 
. ^ , GeV^ range, and individual quark and gluon trial distributions given as functions of x. The distributions are evolved 
■ to larger using the coupled integral-differential DGLAP equations [H H, and the results used to predict the 
5^ I measured quantities. The final distributions are then determined by adjusting the parameters in the input to obtain 
a best fit to the data. For recent determinations of the gluon and quark distributions, see [1, H, @, 0]- 

This procedure is rather indirect, especially so in the case of the gluon: the gluon distribution G{x,Q^) docs not 
appear in the experimentally accessible quantity i^^^(a;, Q^), and is determined only through the quark distributions 
in conjunction with the evolution equations. It is further not clear without detailed analysis how sensitive 

the results are to the parametrizations of the initial parton distributions, or how well the gluon distribution is actually 
determined. 

In two recent papers [Til . [T^ papers, we presented a new method for determining the gluon distribution function 
G{x, Q^) = xg{x, Q^) in leading order (LO) in QCD directly from a global parametrization of the data on F2^{x, Q^), 
and applied the results to obtain sensitive tests of some existing parton distributions. The method uses only the LO 
DGLAP evolution equation [J, [2, lS] for _F-J^(x, Q^) in the usual approximation in which the active quarks arc treated 
as massless. In contrast to previous methods for determining G{x, Q^), it does not require knowledge of the separate 
quark distributions in the region in which structure function data exist, nor does it require the use of the evolution 
equation for G{x,Q^), both considerable simplifications. 

We generalize here to the case of massive quarks using the same methods and the less-stringent approximation 
that mass effects are neglected only in terms found to be small. The additional complications are then minimal. 
Furthermore, to good approximation, detailed knowledge of the individual quark distributions is still not required. 
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In the following Sections, we present two analytic methods that determine G{x, Q^) directly from Q^). The 

first involves an inhomogeneous second-order linear differential equation for G{x,Q^) derived from the LO DGLAP 
equation for the evolution of ^^^(cc, Q^) as a function of Q^. The inhomogeneous driving term in the equation is 
determined entirely by F2^{x, Q^) and its derivatives, and the equation can be solved explicitly. The second method 
involves a somewhat unusual application of Laplace transforms to the DGLAP equation, and leads to equivalent 
results without the intervening differential equation. 

Our method, in either form, determines the LO gluon distribution function G(x, Q^) in terms of the measured 
structure function F2^{x,Q'^). The result is unique, within experimental uncertainties, in the domain D^XjQ^) in 
which there are experimental data and requires no assumptions about the initial shapes of the gluon or the individual 
quark distributions. Furthermore, it is not necessary to use the DGLAP evolution equation for G{x, Q^) as long as 
one remains in the experimental region, though it provides a useful consistency check. These LO results are, of course, 
uncertain at the level of the next-to-leading (NLO ) contributions to the DGLAP equations; these corrections will be 
considered elsewhere. 

As an application of the method, we construct an accurate new global parametrization of the complete set of ZEUS 
data on F:^^{x,Q'^) and use this to determine the 5-quark gluon distribution, G{x,Q^), for massless u, d, s 

and massive c, b quarks. We then discuss in detail the mass effects evident in a comparison with the results for 
massless quarks, including their interesting dependence on . 

We note that the same approach can be used to relate G{x,Q^) directly to data on the remaining DIS structure 
functions F^^^y -^2(3)' ^2(3)' ^'^'^^2(3) -^'-'^ neutral- and charged-current scattering,. 

It is still necessary to obtain the individual quark distributions in order to predict some quantities, or to check 
the accuracy of our approximations. As noted above, this is traditionally done using simultaneous fits to all quark 
and gluon distributions and the complete set of coupled DGLAP equations, fits which may use jet or other quark- or 
gluon-dependent data. However, our results show very clearly the direct connection of G{x,Q^) to the accurate deep 
inelastic scattering data and should provide a useful check on the gluon distributions obtained using other methods. 
We show, for example, that the disagreement of our LO results for massless u, d, s and c quarks with the gluon 
distribution obtained in the CTEQ6L fits [l^ results from the failure of the latter to reproduce accurately the ZEUS 
results for F^P{x,Q'^). 



II. ANALYTIC TREATMENT OF THE LEADING ORDER GLUON DISTRIBUTION FUNCTION 



A. Treatment of heavy-quark effects 



In the presence of heavy quarks, here taken as c and b quarks, the DGLAP evolution equations for the quark 
distribution functions qi{x, Q^), i G u,u, d, d, s, s, c, c, 6, b, must be modified to take into account the effect of production 
thresholds for pairs of heavy quarks, assuming none are present in the initial hadron, and mass effects in the parton 
splitting functions ; each affects the evolution of the parton distributions. We will treat these effects using 
a simplified version of the method introduced by Aivazis, Collins, Olness, and Tung (ACOT) [l^. In this method 
[5l.ll6l.[l7l.[l8lfl9t. a heavy quark has its usual Bjorken scaling variable x replaced in its splitting functions Kgq. [xj z) 
and Kq.q.{x/z) by the "slow" scaling variable Xi = x(l -I- AMf /Q^) that appears naturally. Further mass effects in 
the splitting functions are ignored, so those functions retain the forms used for massless quarks. The integrations in 
the evolution equations are then taken to run from z = Xi to z = 1 rather than z = x to z = 1, thus imposing the 
threshold condition x < Xi < 1. (There are some additional changes if one is interested specifically in the production 
of heavy quarks near threshold [l9| , but these are not important for inclusive deep inelastic scattering and will be 
ignored here.) 

With the modifications above, the evolution equation for a heavy quark becomes 



91ng2 



47r 



dz r 



q,{z,Q^)Kqq (^) +giz,Q^)Kgq (^) 



Multiplying by the squares of the quark charges e| and summing over the quarks i 
2\ 1^ _ ^ _ efqi{x, Q^) on the left, we find that 



w, u, d, d, s, , 



(1) 



to get 



IdF^^ix^i 



51nQ2 



47r 



-jG{z, 



(f) 



(2) 



All terms in the second sum on the right involve G{z,Q^) — zg{z,Q'^), with only the arguments in the splitting 
functions and the ranges of integration modified in the heavy quark terms. This will cause no problems as will be 
seen below. 
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The heavy-quark terms in the first sum on the right do cause a problem: because of the different arguments in the 
the splitting functions and the different ranges of integration, the sum is not expressible in terms of as would be 
the case for massless quarks. This becomes clearer if we introduce the parameters rji{Cf') = 1 + {AMf /Q"^) and change 
the variable of integration in each term from z to z' = z/rji^ and then drop the prime for simplicity. The lower limit 
of the integral becomes x and the upper limit can be extended from l/rji to 1, since qi{z) = for z > 1. The quark 
term then becomes 



as 
An 



1 dz 



%z,Q^)Kq 



1 dz 

^2 2, shifted 



iz,Q')K,, Q 



(3) 



where Kqq is the common splitting function for massless quarks and Fr 



2%iU',Q^) = E.efz<?.(?7.z,g2). When the 
quarks are taken as massless, rji = 1, and the quark sum is just i^^^(z, Q^)/^- This is no longer true when heavy 
quarks are present, and information on the individual quark distributions is clearly needed to evaluate the sum exactly. 

The methods developed below for determining G{x,Q'^) exactly from Eq. ^ require only that F^^ and F^^y^-^^^^^ be 
known. This is of course the case when the quark distributions are known. If they are not, we can still proceed using 
reasonable approximations. The simplest is to replace ?/; by 1 in Eq. ([3]), but not in the gluon terms in Eq. ([2]) where 
threshold effects are important. The quark sum in Eq. ^ then becomes F2^{z,Q'^), and we obtain the simplified 
evolution equation 



IdFlM^ 



An 



(4) 



This approximation should be good. The factor F2^{z, Q^) in the first integral on the right will be evaluated using the 
experimentally determined structure function so automatically includes the threshold conditions for the activation of 
heavy quarks and the expected suppression of the heavy quark distributions. The actual approximation is the neglect 
of the shift from z to rjiZ for the point at which qi is evaluated. This should not be a large effect for sufficiently 
large and qi{x, Q"^) not too rapidly varying in x. In addition, the contribution of the i^^J^-dependent terms in Eq. ^ 
to our final expression for G{x,Q'^) is considerably smaller at small x than that dependent on 9i^^^(a;,Q^)/91nQ^, 
significantly reducing the importance of the approximation. 

We have checked the validity of the approximation using the CTEQ6.5 Q quark distributions which were derived 
using the simplified ACOT method to treat mass effects as is done here. The errors introduced by replacing the 
proper quark sum in Eq. Q by F^P{x,Q'^) is less than 14%, 9%, and 2.5% of F^^ for = 5, 20 and 100 GeV^, 
respectively, for 0.1 > a; > 10~^. Taking into account the expected suppression of the i^^^-dependent terms in the 
final result relative to those dependent on 9F^^(x,(5^)/91nQ^, we expect associated errors in G{x,Q^) considerably 
smaller. They are in fact less than 2-3% in the final G. 

A more accurate approximation for the input that attempts to incorporate the shift from z to rjcZ and rjbZ in the c- 
and 6-quark distribution functions in Eq. ([3]) is to equate those terms to their maxima of 2/3 and 1/6 of the light-quark 
sea contribution to F2^, then shift the arguments appropriately. In particular, taking into account the shifts of the 
heavy-quark variables in Eq. ([T]), we approximate the c and b terms in Fj'^ as 
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xc{x,Q-') « -^F;jP{rj,x,Q^), 



xbix,Q') « —F;iPJrj,x,Q'). 



(5) 
(6) 



Then, noting that the total sea distribution is F^Jl — ^'i^ght + (8/9)xc-|- {2/9)xb, using the approximations above, and 
solving for i^j^ghti fii^d that 



xc{x, Q^) 



3 

4?7c 



xb{x,Q^) « — 



1 + i^n^^c) + 

377c 07^6 



FZlivbX,Q^). 



(7) 
(8) 



Here T{r]) is a scaling operator with the property that T{i])f{x) = f{rix). As will be seen in detail later, T can also 
be viewed as a translation operator in the variable v = ln(l/a;) that shifts v to v — In 77, and can be approximated for 
not too small by the unit operator. 
The problem at this point is that FJ,^^ is not known directly from data. We can estimate it by starting with at 
small x, where the valence-quark contributions are small relative to the sea, and making a smooth extrapolation to 
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zero for a; — > 1 using a reasonable model for the valence distribution. More simply, we can just replace F^,^ by F2'', 
a good approximation at small x and large Q^. In either case, we can use the approximations in Eqs. ([7]) and ([5]) to 
estimate the effect of the shifts from z to ijcZ and 7]i,z in Eq. If we use F^^, we find that 



377c \r]' 



1 



F^''{vix,Q')-F]''{^j,x,Q 



1 



677c \rib 



1 



Frivtx,Q')-F.]P{^tx,Q') 



(9) 



This approximation is quite good at small x. 

Because of the additional complication in intermediate results, we will use the simpler approximation for the 
evolution equation in Eq. and discuss the accuracy of the result later. 



B. Differential equation for G{x,Q^) 

To obtain a differential equation for G{x,Q^), we follow the procedure used in our earlier work but with the 
modified G-depcndent term in Eq. (|4]) instead of that for massless quarks. Writing out the sum and the splitting 
function explicitly for massless u, d, and s quarks and massive c and b quarks, we obtain the equation 



Here Xi = xiji, rji = 1 + {4Mf / Q'^) , and !F]^{x, Q^) is the sum of the i^2''''"dependent terms in Eq. ([4]), 



F7(x,Q2)inL^ 



ain(Q2) 

(mim _ _ 8 . f ^..(,, g.) (1 + ^) ^1 (lla) 

\ Z X / Z 3j J X \ Z / Z \ 



51nQ2 47r [ 3 dz 

The only change if we treat the quark sum in Eq. ([3]) exactly, or in an approximation such as that in Eq. Q , is the 
replacement of F^^ by F^^^y^^^^^^ in all but the first term on the right-hand side of this equation. We will therefore 
assume that J-J'' is known. 

The three terms on the left hand side of Eq. (fTUl) can be combined with a common splitting function by rescaling 
the variables z in the second and third terms to rjcZ and 77^^, respectively. This gives the modified equation 

z'^ \2> 9t]c 9r]b J \ 2 J as X 

As was observed in [ll| , a three-fold differentiation of Eq. ([T^ with respect to x eliminates the integration on the 
left hand side. Dividing out the factor 4/3 in the first term from the sum of the light-quark charges in Eq. pO)) and 
multiplying by a;^, we obtain an inhomogeneous second-order differential equation for G(x, Q^), 

a^'^- 2x^ + 4^ (^G{x,Q') + ^^^G{xe,Q') + ~G{xb,Q^)^ =^3(0:,^'), (13) 



where 
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The functions G{xi,Q'^) in Eq. ((T3l) vanish identically for Xi > 1, conditions that can be made explicit if so desired 
by introducing step functions 0{1 — Xi) with G{xi, Q^) —> 0{1 — Xi)G{xi, Q^). In order to simplify the appearance of 
our equations, we will generally suppress these step functions. 

The function Q3 on the right hand side of Eq. is defined here using the multiplicative factor 3/4 for three 
massless quarks rather than the factor 9/20 used in our treatment of the problem for the 4 massless u, rf, s and c 
quarks in Ref. More generally, for n effectively massless quarks, we define Gn with the factor 3/4 replaced by 

1/^^ef where the sum runs over these "light" quarks and antiquarks. Thus, the function used in [lll | for 4 quarks 
was g4x,Q^) = (3/5)e3(^,Q'). 

The differential equation in Eq. (jl3[) can be put in simpler form by changing to the new variable v = ln(l/a;). Then 
with G(«, Q2) ^ ^(e-", Q2) and G^iv, Q^) = Gsie-^Q^), 



2 1 



1 1 



^+3^ + 4) ( G{v,Q') + =;^G{v + ~G{v -\nf]b) ] ^g3{v,Q''). 



3 ric 



6 Vb 



(15) 



Solving for the w-dependcnt function on the left hand side using standard methods, we find that the exact solution 
for the boundary conditions G(0, Q^) = {dG/dv){0, Q^) = q is 



G{v, Q2) + l-G{v - Inr/e) + ^-G(« - Inry^) = Hsiv, Q^), 
3 77c QVb 



(16) 



where 7^3(1;, Q^) is the function 



A, - A- 



dv' 



Here A± = k zt iuj are the roots of the factored form of the differential operator in Eq. (|15|) : 
{D - X±){D - Azp), D = d/dv, with k = -3/2 and u = V7/2. 
Finally, converting back to x as the variable, we find 



2 1 



1 1 



3 6 77b 



where 



(17) 

4 = 

(18) 
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dz 


A+ 


- A_ 


J X 


z 
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[' dz 




f'- 




L z 


\X 





z\^+ 

X/ 



z\^- 

X/ 



(19a) 
(19b) 



Two steps remain for us to find a solution for G{x, Q^), first, the evaluation of the integral of Gsiz, Q^) which gives 
"Hsix, Q^), and second, the solution of the resulting mixed expression in x, Xc, and Xb for G{x, Q^)- 



C. Calculation of Ha 

After some manipulation starting with the expression for in Eq. (jllbp . we find that t/3(a;, Q^) can be written 



as 



Q^M') = -l-x'l^, - ' '7' +4:^"!^ Z dz—^iz,Q')ln 



z 



4 as dx^ \x (9 In J dx^ \x Jx ^z ' z ~ x 

Because ^ is determined by experimental data with limited accuracy, the appearance of high derivatives (up 
to the fourth) in the expression above could be regarded with suspicion. However, this apparent problem largely 
disappears when we substitute into Eq. (|19ap . We can then integrate by parts to eliminate these derivatives 
as much as possible. The boundary conditions i^^^(l,(5^) — dF2^{l,Q^)/dx = at x = 1 eliminate the leading 
endpoint terms in the integrations. The use of the expression in Eq. (|19aP rather than that in Eq. (|19bp also makes 
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it easy to use the condition that determines the roots in the factored form of the differential operator in Eq. (I15p . 
3A-I- +4 = 0, to eliminate high powers of X± and reduce the results to simpler form. 
Following this procedure, and using the identity for the derivatives of the logarithmic terms given in Eq. (11) of 
Ref. ini, we find that 



(^^'Q )-7— ^7J3nr7^(^'Q ) 



4 as (9 In Q2 



dz ( z 



4 tts dxd In Q"^ 

fc + 3 . z z 
sm ti> in V cos cj in — 

LO X X 

d 



4^ dF^P 
as d\nQ'^ 



-5F^P{x,Q^) + 3x^F^P{x,Q^) 



-12 



-ix 



oz z — X 

2 ^.^^Q2)^^^^iz^^x\T,--^{xM^) 



+8 



dz^ 
dz ( 
z \x' 
^ dz ( z\^ 
z \x ' 



X dx 

20 + 12fc . z z 
sincj In h 12 coscj In — 

UJ X X 



fc + 3 . z z 
sm In — (- cos w In — 

UJ X X 



1 f)piP 

oz' z' 



-dz'. 



(21) 



While somewhat lengthy, this exact expression has the advantage that it determines the gluon distribution directly in 
terms of -Fj^^ and its first two derivatives. Only one double integral remains, a point that simplifies numerical work 
considerably. This would not be the case if we had simply integrated out the explicit derivatives (d/dx)^ that act on 
jx in Eq. (|14p and had then used the expression for in Eq. (|llap . 

The largest terms in 7^3 at small x are those in the first two lines. The leading terms in Eq. (pij) are local and 
can be evaluated directly at the x and of interest using the global parametrization of F^^(a;,Q2). The remaining 
terms involve at most second derivatives. These derivatives are very well determined by the global parametrization 
of F^^ix^Q"^^. The global nature of the fit is essential here: the derivatives would be considerably less certain if we 
attempted to determine them using only local data. 

The problem now is to evaluate the integrals in Eq. (PT|) . These include the region xp < x < \, xp ^ 0.09, not 



covered in the global fit to the F'^'^ data of Berger, Block, and Tan [20[, updated later in Sec. IIII Al of this paper. 
In the absence of a global fit that includes data for xp < x < 1, the integrals can still be evaluated approximately 
by using a reasonable extension of F^^ over that region which satisfies the boundary conditions for x 1, matches 
smoothly to the existing fit for a; < xp, and is consistent with the existing data. The results in the small- a; region 
X <^xp turn out to be insensitive to the form of the extension. 



'H^(x, Q ) can be converted to the v form "^3(1;, Q ) by the substitution x = e 
the integrals are evaluated. 



ln(l/x) either before or after 



D. Laplace transform method 

Before considering the solution of Eq. p6p , we will rederive it directly using a Laplace transform method given in 
which avoids the intermediate steps involved with the differential equation approach, and again shows that high 
derivatives of the structure functions do not appear in Ti^iv^Q^). 

We will start with Eq. p2)) . Dividing out the sum of the squares of the light quark charges in the first term and 
multiplying the equation by x, we write this as 

r^0(.,Q^)^f^^^i^^) (22) 

where 

0(z, Q2) = G{z, Q') + l-Giv,z, Q2) + l-Gimz, Q^). (23) 

Changing variables from x and z to w = ln(l/a-) and w = ln(l/z), Eq. and Eq. can be rewritten as 

/ Q{w,Q^)h{v-w)dw^ f{v,Q^), (24) 
Jo 
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where we have introduced the notation F{v, Q^) for the u-space form of any function F{x, 



The functions in Eq. ([24|) are then 

h{v) = e 



F{v,Q^)=F{e-\Q^). 



2e 



l{v, Q^) = G{v, + l-G{v - In 77,, Q^) + —Giv - Inijb, Q^), 
3 ?7c 6 T]b 



3 47r 

4 a. 



(25) 

(26a) 
(26b) 

(26c) 



Equation (|24p relates /(?;, Q^), a quantity determined by experiment, to the convolution of g, a function dependent 
on the desired function G{v,Q^), with h{v), which is simply e"" times the gluon LO splitting function for the 
production of a quark from a gluon, expressed in u-space. 

The convolution theorem for Laplace transforms applied to transformable functions p and q states that 



C 



p{w)q{v — w)dw] s 

C~'[C[p-s]C[q-s]-v\^ r p 
JO 



^P; s]C[q;s], 
j{w)q{v — w)dw. 



(27a) 
(27b) 



Applying the first form to Eq. ([M)) . we find that the product of the Laplace transforms of the factors in the 
convolution integral is equal to the Laplace transform of /, or solving for the G-dependent factor, that 



C[q{v,Q^)-s\ = [C 
Thus, inverting the Laplace transform on the left, 

0(z;,Q2)^£-i 



h\ s 



C 













[(^ 


h] s 




/>,Q');s 





(28) 



(29) 



The rightmost factor in this equation is not known analytically, so we cannot invert the Laplace transform analyt- 
ically to determine g{v, Q^). However, the Laplace transform of h is easy to calculate. 



C h; s 

as is the (singular) inverse transformation of (^C 

"(s + l)(s + 2)(s + 3) 



3s + 4 



(s + l)(s + 2)(s + 3) 

-1 



(30) 



h: s 



c- 



s2 + 3s + 4 



3S{v) + (5' (w) - e''^ ( - sin ujv + 2 cos ujv 



(31) 



Here fc, uj are the real and imaginary parts of the roots X± ~ k±iuj of the polynomial s^ + 3s + 4 in the denominator on 
the left hand side of Eq. (1311) . The A's appeared earlier as the roots of the differential operator [{d/dv)^ + 3{{d/dv) + 4] 
in Eq. (fT5|) as noted in 

We can now use the second form of the convolution theorem in Eq. (j27bp to evaluate the right hand side of Eq. 
(P^l as a convolution, and find that 

2{v,Q^) ^ 3f{v,Q^-) + -^f{v,Q^) ~ /(u;,Q2)eM— ) Q. smiu{v ~ w) + 2 coslj{v - w)^ dw. (32) 

That is, using Eq. (|26b)) . 

Giv, Q^) + l-G{v - In 77,, Q2) + \-G{v - Injy^, Q^) = n^{v, Q^), (33) 
3 77c 6 77h 
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where 



4 Qfs av 



g2)e'=("-"'' smuj{v -w)+2 coslu{v - w)^ dw \ . (34) 





This result is equivalent to that in Eq. (|16p . with the expression on the right equal to 7i3(w, Q^)- 
Transforming back to x space, we get 

H3(x,Q2) = ^^|3^7(z,Q2)-x^^7(x,g2) 



V7(z,Q2)g) 



3/2 



3 2^ \ diZ 

sin w In — h 2 cos win— I — ^ . (35) 



OJ X X J z 



The equivalence of this to the expression in Eq. (|2ip may be shown making repeated partial integrations to eliminate 
as many as possible of the derivatives and double integrals that appear through ^-"2^. 

The importance of the Laplace construction is in its avoidance of the intermediate differential equation. Further, 
it can also be applied, usin g th e same formalism, to the DGLAP evolution equation for G{x,Q^) and other structure 



functions, as was shown in [1 



E. Solution for G(a:,Q^) 

In the case of four massless quarks considered in Refs. [ll1.[l^ and early treatments of parton distributions, rjc ~ 1, 
the second term in Eq. is just (2/3)G(a;, Q^), the b quark term is absent, and Eqs. (|18p and (|19bp give the exact 
Uf =4 solution for G{x,Q^) in terms of i^^''(a;, Q^), i.e., 

Giix, Q2) ^ ^-Hsix, Q2) ^ mix, Q2) (36) 
5 



The only further calculation necessary is the evaluation of 7^3 using the global fit to F2^{x, Q"^) — see footnote [25 1. 

The situation is more complicated for massive c and h quarks, and we need to solve Eq. p3p (or equivalently, one 
of Eqs. (flBl) or (flS]) ') for G. We will initially use the equation in the v form, in this case explicitly introducing the 
threshold step functions, i.e., 

G(w, Q')0(w) + \—e{v - \n7^c)G{v - In 77c, O') + \-d{v - \nr^b)G{v - Inr/b, Q^) = B{v)nz{v, Q^). (37) 
3 77c 6 77f, 

The second and third terms on the left arc just translates of the first in u-space. The allowable ranges of the 
translations are limited by the threshold condition; the first argument of G must be be positive, with G{w,Q^) = 
for w < 0. Similarly, T-L^(w,Q^) = for w <{). Having shown these ranges explicitly by introducing the appropriate 
step functions 9{w) in the various terms in Eq. (|37p . we will now suppress them in order to keep our expressions 
simple. 

Let a = (2/377c), {3 = (1/6776), and introduce the translation operator T with the property T(u)f{v) = f{v + u). 
Then Eq. ([57)1 can be written as 

[1 + aT{^\nrjc) + PTi-\nrjb)]Giv,Q^) ^n3iv,Q^), (38) 

so 

G{v,Q^) = [l + ar(-ln77c)+/3r(-ln77b)]"'7i3(w,Q'). (39) 
The translations can be implemented explicitly with our Laplace transform technique using the identity 

C[9{v - w)G{v - w, Q^);s] = e~'"''£[G{v, Q^); s]. (40) 
Taking the Laplace transform of Eq. ([37]) , wc find that the equation can be written in Laplace space as 

1 + Ae-'n-'e + ^e-^""'^) C[Giv,Q^);s]=C[nsiv,Q');s], (41) 

3?7c 0776 / 
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so since we earlier set a = 2/3r]c and (3 = l/Grji,, 



s . 



(42) 



The results in Eqs. ([55)) and are exact — see footnote [2q |. 

To make use of these results, we expand the operators on the right in powers of a and /3, and either use the formal 
properties of T in the case of Eq. ([39]) . or calculate the inverse Laplace transform of the series using the result in Eq. 
(jiO]) in the case of Eq. (j42|), to obtain 



N 



G{v,Q^) - 7i3(i',Q') + ^(-ir[aT(-ln7?,) + /3T(-lnr7b)]"H3(«,g') 



N 



[ae 



TV 



k=0 



= Hsiv, Q2) + ^(-1)" E a'^-'^p'^Hsiv -in- k) Inrj, - k In jy^, Q^), 



(43a) 
(43b) 
(43c) 



where (^) is a binomial coefficient. 

As a consequence of the threshold condition that G{w,Q^) = for the argument w < 0, the summations in Eq. 
(|43cp are finite. Since rjc < rji, for fixed Q^, the sum on A; terminates for fixed v and n at the smallest fc for which 
nlnrjc + k{lnr]f) — In jyc) > v. Similarly, the sum on n terminates at N such that (N + 1) In rye > v. The result is 
equivalent to an iterative solution of Eq. pT]) starting with Go = 7^3. 

Converting back to the more familiar variable x, we have 



G{x, Q2) + l-G{v,x, Q2) + 1-G{i]bx, Q2) = Hsix, Q^), 



with the exact solution 



TV 



G(x, Q2) = iy3(.^, q2) + ^(_1)" ^ c."-fc/3'=H3(?7rS,^a;, Q2) 



A;=0 



(44) 



(45) 



The sums terminate when xrf^{rii,/r]c)^ > 1 (fc) and a;?/" > 1 [n). 

The expression in either Eq. (|43c|) or Eq. (|45|) is exact, but may involve a large number of terms before the series 
cut off exactly for large and x small. For example, for ~ h GeV^, rjc = 2.25, and 14 terms are required for 
the c series to terminate exactly for x = 10~^, but only 5 terms for x — 10~^. For = 100 GcV^, rjc = 1.0625, 
Inr/c = 0.0606, and the numbers of terms necessary increase to 152 and 75. Fortunately, the factors a"~*''/3*'' decrease 
rapidly, and very good accuracy can be attained in the sums for much smaller values of N . 

A different approach that is better for much of the region of interest is to rewrite the operator in Eq. (|39p as 



1 



l + aTc + m l + a + fi l + a'{Tc-l)+ l3'{Tb-iy 
where Tc = T{— ln?7c), ?b — T{— ln?7b), and 

a' = a/(l + a + /3), /?' = /3/(l + a + /3), 



(46) 



(47) 



and then expand the last factor. For Inrj small, T(— Inry) approaches the unit operator, and one obtains a rapidly 
convergent series. Specifically, 



G(«,Q2) = 



1 



TV 



l + a + /3 
1 

1 + a + P 



Mv, Q^) + Y.{-lT[a'{Tc - 1) + 13' [n - 1)]"H3(«, Q^) 

\ n=l ) 

n^iv, Q') + ^(-1)" - W-'Win i^nsiv, < 



n=l 



k=0 



(48a) 
(48b) 
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Alternatively, 

G{x, Q') = (^3(x, Q^) + E(-l)" E Q [a'(r, - l)]"-^-[/3'(T, - l)tn,{x, Q^)^ , (49) 

where T™Ti^{x, Q^) = 7^3(77™^, Q^). The results in these equations are still exact. 

It is easily seen that (T — 1) acts essentially as a derivative operator around a shifted point: 

[T{u) - !)]/(«) = (ri/^H - T-^'\u)) T^'\u)f{v) - /(«' + ul2) - f{v' - ujl) « ufiv') (50) 

where v' = v \- u/2. More generally, [T(m) - 1]" f{v) « u'"/(™)(?; + (m/2)w), where = d"^ f (v) / dv"" . 

In the present case, the function 'Hj,{v,Q^) on which T operates is well parametrized for v > 0.9 (x < 0.1) by 
a low-order polynomial in v [ll[, and (T — l)™7i3(ii, Q^) decreases rapidly with increasing m. The convergence is 
further enhanced at large by the smallness of Inyy^ « {AMf /Q"^). We have found in practice that expansion through 
(T — l)'^ is generally enough to obtain the accuracy needed. 

We note finally that for > AM^ w 20 GcV^, a w 2/3, (3 w 1/6, and 7t:3(v, Q^)/(1 + a + /3) w 7t:5(w, Q^), and 
the expansion starts in the limit of five masslcss quarks, as expected. 



III. APPLICATIONS 
A. Global parametrization of {x,Q^) using ZEUS structure function data 

As an application of the procedures above, we will numerically investigate the effects of the c and b masses on 
G{x, Q^) using an updated version of the global parametrization of the ZEUS data for the proton structure function 
F2^{x, Q^) [l3,[l3| made by Berger, Block and Tan Those authors showed that the data for x < xp ~ 0.09 could 
be parameterized very well as a function of and x with the expression — see footnote (27j , 



F7(x,Q2)==(1-x-) 



Fp 



1 — xp 



A(g2)ln 



Xp 1 



X 1 — Xp 



B(g2)ln2 



Xp 1 — X 
X 1 — Xp 



(51) 



Here xp specifies the location in x of an approximate fixed point observed in the data where curves of F2 {x, Q ) for 



different Q'^ cross. At that point, dF^P{xp,Q^)/d\nQ^ 



for aU g^; Fp 



F2^{xp,Q'^) is the common value of 



The g dependence of F2 (x, Q ) is given in those fits by 



A{Q'^) = flo + ailng2 + a2ln^g^ 
-B(g') = bo + hi lng2+62 In' g'. 



(52) 



The original fits to DIS data [H, [l3| given by Berger, Block, and Tan included data at 24 values of Q^, 
g2 = 0.11, 0.25, 0.65, 2.7, 3.5, 4.5, 6.5, 8.5, 10, 12, 15, 18, 22, 27, 35, 45, 70, 90, 120, 200, 250, 450, 800, and 1200 
GeV', and all x < xp, with the scaling point values xp = 0.09 and Fp = 0.41 assumed to be fixed. The fits were 
performed using the "sieve" algorithm [2ll |. which minimizes the squared Lorentzian 



N 



Klicc- a^) EE ^ In {1 + 0.179Ax?(x,; a)} , 



(53) 



where x^(Q:;a;) = X^il^i ^X?(2;j; a), ^xli^i'iOc) = {[yi{xi]oi.) — yi{xi)] /oi)^ , a is the parameter space vector, and 
yi{xi;a) is the theoretical value of the measured yi at Xi, with measurement error ai, using a '-^Ximax '^^^ 6' 
exclude outliers. 

In this paper, we extend the earlier calculation [2^, again using the functional forms given in Eqs. ((5T|) and (|52p . 
but modifying the fit in two important aspects: First, we fit more ZEUS data — 29 different values — now covering 
a larger virtuality range 0.11 < Q'^ < 2000 GeV^ using all of the published ZEUS data sets [H, [13 with x < 0.09, 
i.e., fitting data for Q'^ = 2.7, 3.5, 4.5, 6.5, 8.5, 10, 12, 15, 18, 22, 27, 35, 45, 60, 70, 90, 120, (150), 200, 250, (350), 
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450, (650), 800, 1200 , (1500) and (2000) GcV^. The additional data sets that are now included in the fit are given 
in parentheses. 

Second, we now fit the scaling point (xp, Fp), leaving the values of and Fp as free parameters, to be determined 
by the sieve algorithm. We now find the best value for the scaling points to be xp = 0.0494±0.0039, Fp ~ 0.503±0.012, 
values significantly different from those assumed in the earlier work [20| . 

The new data set has a total of 210 datum points; using the sieve algorithm [2l| eliminated 6 points whose total 
contribution was 61.0. 

The results of the fit are shown in Fig. [1] where we plot F^^(a;,Q^) vs. x. In order not to visually clutter 
Fig. [1] only 13 of the 29 sets used in the fit are shown, in the range from 0.11 to 1200 GeV^, and Bjorken- 
X range, 10~^ < x < xp = 0.049. We see that the compact fit of Eq. (|5ip gives a good parametrization of the 
ZEUS experimental data for the proton structure function F^^(x, Q^) over the available range of data for small x. It 
should be stressed that all curves, independent of their virtuality, go through the scaling point, the intersection of the 
horizontal and vertical straight lines, providing a powerful constraint on the fit. The fit satisfies this constraint with 
a quite satisfactory goodness-of-fit probability, as seen from Table [H 




FIG. 1: Plots of the fitted proton structure function, F2^{x,Q^) vs. Bjorken x, for virtualities — 0.11, 0.25, 0.65, 3.5, 4.5, 
6.5, 10, 15, 22, 35, 70, 250 and 1200 GeV^ The data are from the ZEUS collaboration The curves are from the fit to 

the full data sample. The vertical and horizontal lines intersect at the scaling point xp = 0.049 and Fp — 0.50 determined by 
the fit to the data. 

The values of the 8 fit parameters, along with their statistical errors, are given in Table HI Also shown is the 
corrected pcr degree of freedom, 1.11, for 194 degrees of freedom. This yields a goodness-of-fit probability of 0.15, 
a reasonable value for this much data. 

Our fit to the data on F2^{x, Q^) is so far restricted to the region x < xp; we have not attempted to fit the (less 
accurate) DIS data for x > xp. However, the integrals over F2''' that appear in our expressions for Tisix, Q^) in Eq. 
(|2ip . or in its Laplace equivalent in Eq. (|33p . extend up to a; = 1, so include the interval Xp < x < 1. In the absence 
of a global fit to the data in this region, we will simply extend the parametrization, piecewise, to the large- a: region. 
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TABLE L Results of a 8-parameter fit to F2^{x, Q^) structure function data [l^. using the x and behaviors of Eq. (jSlf) 
and Eg. (|52|) . with in GeV^. The renormalized Xmin per degree of freedom, taking into account the effects of the Axi^ax ~ 6 
cut [21(, is given in the row labeled 7^ x Xmin/d.f. The errors in the fitted parameters are multiplied by the appropriate r^-:>[2l|. 



Parameters Values 





-7.828 X 10"^ ±5.19 x 10" 


-3 


dl 


2.248 X 10"^ ± 1.47 x 10" 


-3 


(12 


2.301 X 10"" ±4.88 x 10" 


-4 


bo 


1.313 X 10"^ ±6.99 X 10" 


-4 


bi 


4.736 X 10"^ ±2.98 x 10" 


-4 


&2 


1.064 X 10~^ ±3.88 X 10" 


-5 


xp 


0.0494 ± 0.0039 




Fp 


0.503 ±0.012 




Xmin 


193.19 




1^ X Xmin 


215.3 




d.f. 


194 




Tl X xiin/d.f. 


1.11 



using the form 



where the exponent n{Q^) is determined by requiring that the functions in Eqs. (jSip and (j54p and their first derivatives 
with respect to x and match at a; = xp. The results are reasonable, and give the required parametrization of 
F2^{x,Q'^) over all x. The results for G{x,Q^) in the small- a; region x <^ xp are insensitive to the contributions to 
the integrals from the interval xp < x < 1, hence, to the details of the extension of the parametrization. 



1-x 



B. Evaluation of Hsix, Q^). 

We now have a complete parametrization of F2^{x, Q^) using Eq. ([51]) and Eq. (|52p . for small x, and Eq. ([54]) for 
large x, so we can evaluate Tl^{x, Q"^), using Eq. (|llap or (|llbp . In these calculations, we use the LO form of as{Q^), 

with as(M|) = 0.118 and A5 = 87.8 MeV. This is matched to the expression for n/ = 4 at ^ 7\/2 ^j^j^ jyj^ ^ 4 5 
GcV, giving A4 = 120.4 MeV, which is matched in turn to the expression for nf ^ 3 a.t = with Mc = 1.3 GcV, 
giving A3 = 143.6 MeV. Finally, wc find 7^3(2;, Q^) by inserting Q^) into Eq. (gT]) or Eq. and evaluating 

the integrals numerically. 

We find that we can fit the resulting T-L^{x, Q^) very well for small x with an expression quadratic in both In(Q^) 
and ln(l/a:) , with 

n3ix,Q^) = -2.94- 0.359 ln(Q2)-0.1011n^(Q2) 

+ (0.594- 0.0792 ln(Q2) _ 0.000578 ln2(Q2)) ln(l/a;) 

+ (0.168 ±0.1381n(g2)+ 0.0169 In2(g2))ln2(l/.T), 0<x<xg, (56) 



where xq = 0.06. For the extension to large x, we use 



■H3{x,Q') ^ GixG,Q')( — ] (^-^] : XG<x<l, (57) 

\xgJ \1-xgJ 

where we determine p{Q^) by matching the first derivatives of Eq. ([55]) and Eq. ((57|) a.t x = xg- 
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The function 7^3(2;, Q^) is determined by the measured Q^), so is taken as fixed when we determine the 

gluon distribution G{x,Q^). The transformation to G{x,Q'^) depends on nj and the quark masses. However, in the 
simple case of massless quarks, 'Hz{x,Q'^) is just the gluon distribution Gz{x,Q'^) for n/ ~ 3 massless quarks. The 
massless distributions for Uf = 4 and Uf = 5 massless quarks are just 3/5 and 6/11 of H3{x,Q^). We expect the 
n/ = 3 and Uf = 5 massless distributions to represent the extreme cases of G{x,Q^) for massive quarks, with the 
massive gluon distributions approaching G3{x,Q^) from below as 0, and G5{x,Q^) from above as 00. 

This will be evident later in Fig. [S] 

We now turn to applications of our results. 

80 
60 

a 40 
20 



10"^ lO"'* 0.001 0.01 0.1 

X 




FIG. 2: We compare LO GcTEQ6L(a;, Q^) (the thin curves ) with our massless LO G{x,Q^), for = 5 GeV'^ and Uf — 4 
(red solid curve); — 20 GeV^ and Uf — A (green dashed curves); and at = 100 GeV^ and nj —5 (blue dotted curves). 
The G4 plots for = 5 and 20 GeV^ were made using G4 = (3/5)7^3, and the G5 plot for = 100 GeV^ was made using 
G5 ~ (6/11)7^3, with 7^3 determined from the ZEUS structure function data as described in the text. 



1. Uncertainties in the LO massless gluon distributions 

We immediately deduce from Eq. ((35)) that the fractional uncertainty in Ti.3{x, Q^) is determined by the fractional 
statistical error in J^2^{x,Q'^), where JF^^(x,(5^) is given by Eq. (jllap . As observed earlier in Section FlI Al in our 
discussion of Eq. the contribution to J^2^{x,Q'^) is dominated by the term dF2^{x,Q^)/dlnQ^ . Since the LO 
gluon distribution Gn(x, Q^) for n massless quarks is a constant multiple of Tlsix, Q^), we can estimate the fractional 
statistical error in G„ simply as 

AG'.(x,Q^) ^ An3{x,Q^) ^ AF^Pjx.Q^) _ 

G„(x,Q2) n3{x,Q^) ^ F^Pix,Q^) ' " ^^''> 

with AGi and AF2^ the errors in the functions G„ and the structure function i^^^, respectively, with the errors in 
the structure function being the statistical errors of the fit parameters. Using standard error analysis techniques, we 
calculated the fractional statistical errors AF2^{x,Q'^)/F2^{x,Q'^) from the fit parameters shown in Table HI using 
the correlations (not shown) as well as the diagonal elements of the mass squared matrix of the fit. For x < 0.01, they 
were found to be: < 2% at = 5 QeV^, < 2.5% at = 20 GeV^, < 2.7% for = 100 GeV^. We note that 
for LO massless gluons, the error estimations are straightforward; the statistical error is effectively the total error for 
low X — no significant approximations arc made in this regime 

For comparison, the CTEQ6M total gluon fractional uncertainties Q for = 10 GcV^ range from ^ 20% at 
X = 0.0001 to ^ 10% at a; = 0.01, errors that are considerably larger. Presumably, this is because of the indirect 
method that CTEQ employed to get their gluon distributions, solving the coupled DGLAP evolution equations 
simultaneously for the individual quark and gluon distributiions, and to their inclusion of other data incompatible 
with the ZEUS data. Their techniques for estimating errors in this more complicated situation arc also different. 
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2. Comparison of the massless nj — 4 and Uf — 5 LO gluon distributions to CTEQ6L 



The LO CTEQ6L gluon distributions were derived using the complete set of quark and gluon evolution equations, 
and a scheme for taking mass effects into account for the c and b quarks that depends only on Q^, and does not mix 
the behavior in x and as happens in the simplified ACOT scheme. The CTEQ6L treatment of ^^(Q^) is the same 
as used here, and as discussed in [iJl, the scheme used for treating masses is such that we can directly compare our 
calculation of G{x, Q^) for massless quarks using ny = 4 to the CTEQ6L result in the region < < M^, and 
our Uf = 5 massless calculation to GcTEQeL for > M^. The results should agree. We show this comparison in 
Fig. pi for = 5, 20 and 100 GeV^. The thin curves are those for LO GcteQ6L [IBIj taken from the Durham web 
site [22|, and the thick curves are for our massless solutions. The solid red curves are for = 5 GeV^, the dashed 
green curves for = 20 GeV^, and the dotted blue curves for = 100 GeV^. 

These plots of G should be identical over the region of x and covered by the experimental ZEUS data 
provided we both fit those data. The minimum x-value experimentally observed for a given are a;min(5 GeV^) = 
1 X lO"'', a;min(20 GeV^) = 5 x 10"'^, and Xmin(100 GeV^) = 2 x 10"^. We sec from Fig[2]that, for each of the values 
of plotted, there are significant numerical differences for Bjorken-x greater than .Tmin(Q^), i.e., in the regions in 
which the two determinations of G{x, Q^) should be the same numerically. 

We found earlier [l^l that the CTEQ6L gluon distributions GcTEQ6L(a;, Q^) were consistent with the structure 
functions F2^{x,Q^) calculated from the CTEQ6L quark distributions using Eg. ([55]) . We thus conclude that the 
CTEQ6L quark distributions must not give a very good fit to the ZEUS proton structure function F2^{x, Q^), which 
was the starting point of our calculation. This is indeed the case, as is seen in Fig. [31 where thir results systematically 
miss the experimental data points. We note that the CTEQ6L fits also use other data, suggesting that incompatibilities 
in the distributions required by different data sets may affect their fit to 




FIG. 3: A plot of the proton structure function F2^{x,Q'^) vs. x, for several virtualities. We compare our proton structure 
function fit — the thick curves, using the results of TableUin Eq. I|51ll and Eq. (|52p — to the {x,Q^)cTEQeh, the thin curves, 
and to the ZEUS experimental proton structure function data [3, [l3]. The solid (red) curves are for Q2=4.5 GeV^ the 
dot-dashed (blue) curves are for GeV^ and the dashed (green) curves are for Q^=90 GeV^. The intersection of the 

horizontal and vertical lines indicate the scaling point, shown in Table HI The CTEQ6L values were constructed from their 
quark distributions [2^. 

To study this disagreement in more detail, we constructed i^^^(a;, Qf)cTEQ6L = J2i ^i^1i,CTEQ6L{x, Q^) , using the 
CTEQ6L [131 quark distributions taken from the Durham web site [221 ■ In Fig. |3lwe show plots of F2^{x, Q^) vs. x 
for our fit — using the results of TableUin Eq. ([5T|) and Eq. (|5^ — as the thick curves and the CTEQ6L fits as the thin 
curves, comparing them to the experimental ZEUS proton structure function data [isl. [l^. The solid red curves are 
for = 4.5 GeV^, the dot-dashed blue curves are for — 22 GeV^, and the dashed green curves are for = 90 
GeV^, allowing us to compare virtualities that are very close to those used in Fig. |2l Clearly, the CTEQ6L curves 
are significantly lower than the experimental data, as well as having a different dependence, resulting in a poor 
fit to the experimental ZEUS data. Conversely, our fit is a good representation of the proton structure function data. 
This difference explains the discrepancy, in the appropriate x regions, between the numerical values of the different 
gluon distributions shown in Fig. |2l 

We re-emphasize that our method connects G{x,Q^) in LO directly to the structure function data, without the 
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need to involve the individual quarks distributions. Other structure function data can be used the same way to derive 
independent expressions for G [l^l- Since the evolution of G{x,Q^) is built implicitly into the data, we do not need 
to use the gluon evolution equation to obtain our results. However, we could still use the gluon evolution equation to 
check the consistency of our results, e.g., through a determination of Fs [l^l, giving a test of the adequacy of the LO 
approximation of the DGLAP equations, a work that is in progress. 

3. Mass effects in nj =4 distributions 

To investigate the effect of masses on G{x, Q^) in the simplified ACOT scheme [l^, we have done model calculations 
for = 5 and 20 GeV^, including only the c quark, and using the approximation of replacing -F'2'^sijiftod(2^i Q^) Eq. 
([3]) by F2^{x, Q^) as in Eq. The results of the calculations are shown in Fig. [H 

To check the approximation, we have repeated the calculation with i^^^ replaced by F2^{x,Q^) + AF2^{x,Q^), 
with the shift correction AF2^ evaluated using the CTEQ6.5 0] quark distributions. The results of the modified 
calculations cannot be distinguished from those in Fig. [4] on the scale of the figure, and are not shown. 

The striking mass effects evident in the differences between the solid red (massive) and dashed blue (massless 
Uf = 4) curves in the figure result entirely from the mass dependence on the left-hand side of Eq. ([44]) , with the b 
term absent. The convergence of the (T^, — 1) expansion for G{x,Q^) in Eq. (j48ap is very rapid, and the magnitude 
of the difference between the massless and massive curves results largely from the overall normalization factors, 3/5 
and 1/[1 + (2/3?7c) in the equations which relate G4 = H4 and G to TI3, Eqs. ([55]) and 

C. Evaluation of the Uf — 5 LO gluon distribution G{x,Q^) for massive c and b quarks 

After evaluating the LO Uf — 5 gluon distribution G{x,Q^) numerically for massless u, d, s and massive c and b 
quarks using the methods described in Sec. Ill El we found that we could obtain an excellent fit to its x dependence 
for small x and fixed using a quadratic expression in In 1/x or v, just as in the case of TI3 and the related massless 
distributions. However, the shape of G{x, Q^) as a function of is far from quadratic, and it requires a much more 
complicated power series in \n{Q^) to fit: 

G{x, Q^) = -2.65 - 0.3671n(g2) + 0.146 In^iQ"^) ~ 0.0321 In^iQ"^) + 0.00160 In^(Q^) 

+ (0.584 + 0.0416 ln(Q2) _ 0.O666 In^iQ"^) + 0.0101 In^iQ^) ~ 0.000490 In^(Q^)) ln(l/x) 
+ (0.155 + 0.123 ln(02) _ 0.0121 ln^(Q2) + 0.00270 ln^(Q2) - 0.000111 ln^(Q2)) ln2(i/a;) 

for < X < = 0.05. (59) 

The x-dependence of Eq. (|45p . the gluon distribution, G{x,Q^) for n/ = 5, with massless u,d,s and massive c,b 
quarks is plotted in Fig. [5] for representative values of Q^, namely = 5, 20 and 100 GeV^. The red dashed curve 
was obtained for = 5 GeV^, the green dot-dashed curve for = 20 GeV^ and the blue dotted curve for = 100 
GeV^. Again, as was true for the massless cases, we see that they have a quadratic dependence on ln(l/a;), with 
the actual value of G{x,Q^), at any x, lying between the massless cases for nj ~ 3 and jif = 5. For x < 0.01, the 
systematic error estimates from the approximations that were actually used (see Section [II A[) are comparable with the 
statistical errors discussed in Section [III B 11 yielding total errors in the 3 — 4% range. If we had employed the more 
accurate approximations discussed in Section III Al the errors would have effectively been reduced to the statistical 
errors of the fit to the ZEUS data [H, [H for F:^P{x, Q^). 

It is evident from the monotonic decrease of G{x, Q^) with increasing x, or, equivalently, the decrease of G{v, Q^) 
with decreasing v, that the values of G{r}x, Q^) or G{v — Inry) at the shifted points that appear in the basic equations, 
Eqs. and ([T5|) are less than G{x, Q^) = G{v, Q^). As a result, using the fact that i]i > 1, 

Giv,Q^) + l-Giv~\nrjc) + l-G{v^liir]b) < (ll/6)G(f , Q^). (60) 
3 7]c 6 r]b 

The function on the left-hand side of this inequality is equal to the fixed function 7i3{v, Q^), so wc find that 

G(«, Q2) > (6/11)71:3(1;, Q') = G^{v, Q"). (61) 

Also, since 

G(t;, g") + I -G(z; - In ryj + \-G{y - In 775) > G(7;, Q^) , (62) 
3 ?7c 6 ?76 
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FIG. 4: Plots of G{x,Q^) vs. x for u, d, s, and c quarks, with the c quark treated as massless (blue dashed curves) and as 
massive (red solid curves): (a) = 5 GeV^; (b) = 20 GeV^. The input function H3{x,Q^) in Eqs. and based 
on the new global fit to F2^{x,Q'^) is given in Eq. (|56p . The changes in G that result from using the function -F2'shiftcd(^' Q^) 
with a shifted argument for the c contribution are small, and the resulting curves cannot be distinguished from those shown on 
the scale of the figure. The differences between the dashed and solid curves result entirely from the mass dependence on the 
left-hand side of Eq. (|44)) . 

we find that G{v, Q^) < Hsiv, Q^) = Gsiv, Q^). 

Thus, for massless u, d, s quarks and massive c and b quarks, G{v,Q^) = G{x,Q^), must always lie between the 
limiting expressions Gnf{x,Q^) for Uf = 3 and Uf = 5 massless quarks. In particular, G approaches G^{x,Q^) from 
above for Q"^ — + oo. In 77 ~ 4A'/^/(5^ 0, and G3{x,Q^) from below for 0, where rjx 1, G{r]x,Q^) — *■ 0, 

and the b and c terms in Eq. (|44p vanish in succession. The limits are, of course, expected: quarks' masses become 
irrelevant for ^ AAlf, and, conversely, heavy quarks are not excited at fixed x for small enough that the 
threshold condition Xi = rjiX < 1 for pair production from a gluon cannot be satisfied. 

This behavior of G{x,Q'^) for massive c and b quarks is shown in Fig. [5] over a very wide range of for two 
representative values of a;, x = 0.0001 and x = 0.01. For comparison, we have also plotted — for the same x — the 
dependence of the massless n/ = 3 gluon distribution G3(x, Q^) as the thin dot-dashed blue curve, the massless 
Uf = A gluon distribution Gi{x,Q'^) as the thin dotted green curve, and the massless nf = 5 gluon distribution 
G5{x, Q^) as the thin soHd red curve. 

From Figure [51 we see that all of the massless distributions rise smoothly with In as expected for a dominant 
In^ behavior as we found for T-L^{x,Q'^), Eq. ((56)) . However, the thick black dashed curve, the nj = 5 distribution 
G{x,Q'^) for massive c and b quarks is clearly far from quadratic. As 0, G{x,Q^) approaches the 3 massless 

quark distribution from below, differing from G3 by only a few percent at ~ 1 GeV^. At ~ 50 GeV^, G crosses 
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FIG. 5: Plots of the final 5-quark gluon distribution for massive quarks, G{x,Q^ 



We show the LO gluon distributions, 



G{x, Q ) for 71/ = 5 with massless u, d, s and massive c and 6 quarks for; Q = 5 GeV , the dashed red curve; = 20 GeV", 
the dot-dashed green curve; and for = 100 GeV^ , the dotted blue curve. The plots were made from Eq. (|45p which used the 
value of Tl3{x,Q^) obtained from the new parametrization of F2{x,Q^) in Eqs. (|5H) and (|52|l with the parameters in Table |T] 
and the LO (n/ =5) form of as{Q^) from Eq. ([55|l . 



the 4-quark case, Gi. For ~ 1000 GeV^, well over the 6-quark threshold, G approaches the 5-quark massless 
solution G5 from above, all as expected from the discussion above. In between, we have a rather complicated 
dependence. We find these patterns to be true for all small x. 

The deviation of the dependence of G{x,Q^) from the smooth behavior of 7^3(2:, Q^) is largely accounted for 
by the prefactor l/(l + a + /3) in Eq. The series this multiplies converges very rapidly and is dominated by the 

leading term. Recalling that rji = 1 + (4Mf/Q^), we see that the prefactor varies from 1 for « 4M^ to 6/11 for 
» 4A/^, and the leading approximation to G has the limits discussed above. 

The term 1/(1 + a + /3) in fact plays the role of a ratio of the sum ef — 4/3 for the three massless quarks u, d, s 
and their antiquarks used to normalize and its integral TLz, to a sum of squares of charges for all the quarks. 



(1 



8 1 
9^ 



2 1 
9^ 



(63) 



The terms 2e^ appear in this sum with weights 1 /r]i which we may associate with the degree to which quark i is active, 
varying from no excitation for << AMf to complete excitation for >> AMf. G{x, Q^) is then given in leading 
approximation for massive quarks by the same expression as for massless quarks, G{x,Q^) « (4/3)(l/e^fj)H3(a;, Q^). 
To the extent to which we may usefully talk about an effective number of active quarks, «/_eff, the preceding argument 
suggests that it should also be defined in terms of the weights l/rji, with n/^off = 3-1- (l/?7c) + i^/vt)- 



IV. CONCLUSIONS 



We have demonstrated that a parametrization of the ZEUS experimental data [H, [3l on the proton structure 
function F2^{x,Q'^) as a function of x and in the domain Vix^Q^) (shown in Fig. [1]) is all that is needed 
to obtain an analytic solution (in the same domain V) for the LO gluon distribution Gnf{x,Q'^) — xg{x,Q^) for 
Uf massless quarks, since Gnf is a numerical multiple of a function 'H^{x,Q'^) that is completely determined by 
F2^{x,Q'^). Comparison with CTEQ6L gluon distributions in the same domain — where they should agree — show 
significant inconsistencies; however, these are explained by the fact that the structure function F2^(x, Q^) constructed 
from the CTEQ6L quark distributions disagrees markedly with the ZEUS data in the domain 2?. 

The same procedure, again using only experimental structure function data, also gives an excellent approximation 
to G{x,Q^) when the c and b quarks are properly treated as massive, using a simplified ACOT approximation p^ . 
In that case, G{x, Q^) is the solution of the equation 



G(x, Q2) + l-Gixc, Q^) + l-G{xb, Q^) = Hsix, 



(64) 
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FIG. 6: G{x,Q'^) vs. Q^, in GeV^, for (a) x = 10"* and (b) x = 10"^ . The thin dot-dashed blue curve is the n/ = 3 plot, 
Gs, for massless u,d,s quarks; the thin dotted green curve is the nf — 4 plot, G4 for massless u,d,s,c quarks; the thin solid 
red curve is the n/ = 5 plot, G5 for massless u, d, s, c, b quarks. The thick dashed black curve is G, the nf — 5 distribution for 
3 massless quarks u,d,s and 2 massive quarks c,b, from Eq. (|45|) . For the evaluation of TL'i{x,Q'^) that was used in Eq. (|45fl . 
we have used the as(Q^) of Eq. (|55p . It should be noted that asymptotically, as 0, G{x,Q^) Gz{x,Q^) from below, 

but as 00, G(x,Q^) G5{x,Q^) from above, and doesn't appear to significantly overlap G4 anywhere. 

where rji = l + {AMf /Q'^) and G{r].iX, Q^) = for rjiX > 1. This was derived and solved using two different approaches, 
one based on a differential equation derived from the DGLAP evolution equation for F^^, and and a second which 
uses a Laplace transform method to solve that equation directly. 

Taking mass effects into account, the analytic solution G(x,Q^) for 5 quarks — three massless u, d, s quarks and 
the massive c and b quarks — is given in terms of Tisix, Q^) by 

n—lm—Q ^ / \ /c / \ / / 

^?/r™%" < 1. (65) 

As an application of our approach, we have made an accurate paramctrization of F2^{x, Q^) using all available low-x 
ZEUS proton structure function data [H, [l3]j and have evaluated H3{x,Q^) numerically. We found that 7i3(x, Q^) 
and the massless solutions for G for small x are expressible numerically as relatively simple quadratic functions of both 
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ln(l/s) and In(Q^). In the massive quark case, G{x,Q'^) is still quadratic in \nl/x, but is a much more complicated 
function of In(Q^) which is bounded from above as — *■ by the massless 3 quark distribution G3{x, Q^), and from 
below as cx) by the massless 5-quark distribution G^^x, Q^)- The same methods can be used to extract G{x, Q^) 

in LO from data on other nuclcon structure functions in DIS. 

We are currently working on NLO effects on gluon distributions for both massless and massive quarks, as well as 
checking the consistency of the DGLAP evolution equations for F2^{x, Q^) and G{x, Q^) in LO, using proton structure 
function data. 
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